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^ Abstract 

In this paper, wc contribute operator-splitting methods improved by the 
Zassenhaus product for the numerical solution of linear partial differential equa- 
^ tions. We address iterative splitting methods, that can be improved by means 

1—1 of the Zassenhaus product formula, which is a scxjunential splitting scheme. 

^ The coupling of iterative and sequential splitting techniques are discussed and 

^ can be combined with respect to their compuational time. While the itera- 

tive splitting schemes are cheap to compute, the Zassenhaus product formula is 
00 more expensive, based on the commutators but achieves higher order accuracy. 

CO Iterative splitting schemes and also Zassenhaus products are applied in physics 

and physical chemistry are important and are predestinated to their combina- 
tions of each benefits. Here we consider phase models in CFD (computational 
fluid dynamics). We present an underlying analysis for obtaining higher order 
operator-splitting methods based on the Zassenhaus product. 
. • Computational benefits are given with sparse matrices, which arose of spa- 

^ ^ tial discretization of the underlying partial differential equations. While Zassen- 

^> haus formula allows higher accuracy, due to the fact that we obtain higher order 

^ commutators, we combine such an improved initialization process to cheap com- 

^ putable to linear convergent iterative splitting schemes. 

Theoretical discussion about convergence and application examples are dis- 
cussed with CFD problems. 

Keywords. Operator splitting method, Iterative splitting, Zassenhaus product, 
Parabolic differential equations, Convection-DifFusion-Reaction Equations. 
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1 Introduction 

Our motivation to study the operator splitting methods come from models in fluid 
dynamics problems, for example problems in bio-remediation [IJ or radioactive con- 
taminants [1], [3]. 

While standard splitting methods deal with lower order convergence, we propose 
to a combination of iterative splitting methods with embedded Zassenhaus product 
formula. 

Theoretically, we combine fix-point schemes (iterative splitting methods) with se- 
quential splitting schemes (Zassenhaus products), which are connected to the theory 
of Lie groups and Lie algebras. Based on that relation, we can construct higher order 
splitting schemes for an underlying Lie algebra and improve the convergence results 
with cheap iterative schemes. 

Historically, the efficiency of decoupling different physical processes into more sim- 
pler processes, e.g., convection and reaction, helps to accelerate the solver process 
and is discussed in [18]. 

We propose the following ideas: 

• Iterative splitting schemes are based on fix-point schemes, e.g. Waveform 
relaxation, which linearly improve the convergence order. Based on reducing 
integral operators to cheap computable matrices, they can be seen as solver 
methods, see [TlJ . 

• Zassenhaus formula are based on nested commutators, which are main keys 
to derive higher order standard splitting schemes (e.g. Lie- Trotter and Strang 
splitting). They are simple to compute with their nil-potent structure, see 

m- 

In this paper we study the following mathematical equations: The equations are 
coupled with the reaction terms and are presented as follows. 

dtRiUi + V • V Uj = -Xi Ri Ui + Ai-i Ri-i Ui-i (1) 
+l3{-Ui+gi) inf] X (0,r) , 

Uifi[x) = Ui[x,{)) on , (2) 
dtRiSi = — Aj Ri gi + Aj-i gi-i (3) 
+P{-gi + Ui) inO X (0,r) , 

9i,o{x) = gi{x,0) on n , (4) 
i = 1, . . . ,m , 
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where m is the number of equations and i is the index of each component. The 
unknown mobile concentrations Ui = Ui{x, t) are considered in x (0, T) C M" x M"*", 
where n is the spatial dimension. The unknown immobile concentrations gi = gi{x, t) 
are considered in Q x (0,T) C M*^ x M"*", where n is the spatial dimension. The 
retardation factors Ri are constant and Ri > 0. The kinetic part is given by the 
factors Aj. They are constant and Aj > 0. For the initialization of the kinetic part, 
we set Ao = 0. The kinetic part is linear and irreversible, so the successors have only 
one predecessor. The initial conditions are given for each component i as constants 
or linear impulses. For the boundary conditions we have trivial inflow and outflow 
conditions with Uj = at the inflow boundary. The transport part is given by the 
velocity v G M" and is piecewise constant, see and [6j. The exchange between 
the mobile and immobile part is given by /3. 

The outline of the paper is as follows. The splitting- methods are discussed in Section 
[2] In Section [3j we present the numerical experiments and the benefits of the higher 
order splitting methods. Finally, we discuss future work in the area of iterative and 
non-iterative splitting methods. 



2 Operator splitting methods 



We focus our attention on the case of two linear operators (i.e., we consider the 
Cauchy problem): 



= Ac(t) + Bc(t), with t G [0,r], c(0) = cq, (5) 

ot 

whereby the initial function cq is given and A and B are assumed to be bounded 
linear operators in the Banach-space X with ^, i3 : X — t- X. In realistic applications 
the operators correspond to physical operators such as convection and diffusion op- 
erators. We consider the following operator splitting schemes: 



2.1 Iterative Operator Splitting 

Iterative splitting with respect to one operator 

^^=Ac,{t) + Sc,_i(t), with = c^,i = 1,2,..., m (6) 

at 

Theorem 1. Let us consider the abstract Cauchy problem given in Then, we 
the one-side iterative operator splitting method ^ has the following accuracy: 

\\{S,-expi{A + B)T)\\<CT\ (7) 
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where Si is the approximated solution for the i-th iterative step and C is a constant 
that can be chosen uniformly on bounded time intervals. 

Proof. The proof is done for i = 1,2, . . . and with the consistency error the ei(r) = 
c(r) — Ci(r) we have : 



air) = exp(ylr)c(r) (8) 

, s , ds 

IfT, 



+ [ exp(y4(t"+^-s))Bexp(s^)c(r) 
Jt" 

+ / exp(A(r+^ - si))B / exp((r+^ - si - S2)A)B ex.-p{s2A)c{t'^) ds2 dsi 

+ ... 

+ / exp(y4(r+^ - si))B / exp((r+^ - si - S2)A)B exp(s2A)c(t") ds2 dsi + ... 

Jt" Jt"' 

+ / exp(^(r+i - si))B 
Jt" 

• • . / exp((t"+ - > Sj)A)Bcinit{si) dsi dsi-i ... ds2 dsi, 

Jt- ,=1 

c(r) = exp(Ar)c(r) (9) 

+ / exp(A(r+i - s))5exp(sA)c(r) 
Jt" 

+ / exp(A(t"+i - si))S / exp((t"+i - si - sa)^)^ exp(s2A)c(t") dsa dsi 

.7t" Jt" 
+ ... 

+ / exp(^(r+^ - si))B / exp((r+^ - si - S2)A)B ex.p{s2A)c{t'^) ds2 dsi 

Jt" Jt" 

+ ...+ / exp(A(r+i - 
Jt" 

... ' exp((r+^ - Vsj)^)5exp((si(^ + 5))c(t") dsi . . . ds2 dsi, 

Jt" „_i 



We obtain: 

lle^ll < ||exp((A + i?)r)c(t") - 5,(r)Q„,t(r)|| 

< Ct' max II exp((si(^ + B))c{t'^) - Cinit{si)\\ 

sie[o,T] 

< CT'\\eM{r{A + BMt^)-Cinit{r)\\, (10) 
where i is the number of iterative steps. 
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The same idea can be applied to the even iterative scheme and also for alternating 
A and B. 

□ 

Remark 1. The accuracy of the initialisation CmitiT) is important to conserve or 
improve the underlying iterative splitting scheme. 

Here we have the following initialization schemes: 

air) = exp(AT)c(r) ^ llcill < Cr*c(r), (11) 
air) = exp(.4r)exp(Sr)c(t") ^ ||ei|| < Cr^+^^t"). (12) 

2.2 Zassenhaus formula (Sequential Splitting) 

The Zassenhaus formula is an extension to the exponential splitting schemes and is 
given as: 

exp((^ + B)t) = nj^^ exp(ai^t) exp(6iSf)7r^^^. exp(Cfct^) + 0(i'"+i). (13) 
where Cj is a function of Lie brackets of A and B. 

Theorem 2. The initial value problem ^ is solved by classical exponential splitting 
schemes. 

Then we can embed the Zassenhaus formula and improve the classical splitting 
schemes 

exp{{A + B)t) = tt{^^ expiuiAt) exp{biBt)n^^j exp(Cfct'^) + 0(t™+i). (14) 
where Cj is a function of the Lie brackets of A and B. 

Proof. 1.) Lie- Trotter splitting: 

For the Lie- Trotter splitting there exists coefficients with respect to the extension: 
exp((^ + B)t) = exp(At) exp(5t)n^2 exp(C7fct'=), (15) 
where the coefficients Ck are given in |10j . 

Based on an existing Baker-Campbell-Hausdorff (BCH) formula of the Lie- Trotter 
splitting one can apply the Zassenhaus formula. 

2.) Strang Splitting: 

A existing BCH formula is given as: 

exp(^t/2) exp(5t) exp(^t/2) = exp(t5'i + t^S's + t^S^ + ...), (16) 
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where the coefficients Si are given as in |12| . 

There exists a Zassenhaus formula based on the BCH formula: 



exp((A/2 + B/2)t) = exp{Ckt^) exp(^/2t) exp(5/2t), (17) 

and 

exp{{B/2 + A/2)t) = exp{B/2t) exp(^/2t)n^2 exp(Cfct''), (18) 
then there exists a new product: 

nr=3 expiDkt'^) = nr=2 exp(C'fct'^)nr=2 expiCkt''), (19) 
with one order higher, see also [22\ . 

□ 

Remark 2. In the following, we concentrate on the Lie- Trotter splitting with the 
embedded Zassenhaus formula, given as: 

Ezassen,Comp,i{t) = exp{At) exp{Bt), (20) 
Ezassen,Comp,j{t) = exp{Cjt^),for j £2...,i, (21) 

where the sequential Zassenhaus operator 

Ezassen,i(,i") — ^j=lE/Zassen,Comp,j(j') , (2^) 

is of accuracy 0(t*) and i are the number of Zassenhaus components. 

2.3 Embedding Zassenhaus expansion to Iterative Splitting schemes 

In the following we discuss the embedding of the Zassenhaus formula into the iter- 
ative operator splitting schemes. 

Theorem 3. We solve the initial value problem We assume bounded and con- 
stant operators A, B. 

The initialization process is done with the Zassenhaus formula: 

Ci{t) = Ezassen,i{t)co. (23) 



where Ezassen,i{t) is given in (22). 



Further the improved solutions are embedded to the iterative splitting schemes ^ 
and we have after j iterative steps the following result: 

Ci+j{t) = Eiter,j{t)Ez assen, ^{t)co. (24) 

where we can improve the error of the iterative scheme to 0{t'^^^). 
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Proof. The solution of the iterative sphtting scheme ([6]) is given as: 

Ci+j{t) = Eiter,i{t)Cinit,j- (25) 

where Si{t) = Eiter,i{t). 

The initiahzation is given with the Zassenhaus formula as: 

Cinit,j{t) = Ezassen,j{t)co. (26) 

combining both splitting schemes we have the local error: 

ei+j{t) = \\c{t) - Ci+j{t)\\ = \\c{t) - Eiter,i{t)Ezassen,j{t)co\\ < O{f+^)co. (27) 

□ 



3 Numerical Examples 

In this section, we discuss examples to the usage of the embedded Zassenhaus prod- 
uct methods to the iterative splitting schemes. In the first examples, we demon- 
strate somewhat artificially how the proposed Zassenhaus splitting method avoids 
the splitting error up to a certain order. The next examples show the solution of 
partial differential equation which can be improved by the Zassenhaus products. 

In the following, we deal with numerical example to verify the theoretical results. 



3.1 First Example: Matrix problem 

For another example, consider the matrix equation. 



u'{t) 



1 2 
3 



u, u(0) = Uq 



(28) 



the exact solution is 



u{t) = 2(e3* - e-"*)/5. 



-2t\ 



(29) 



We split the matrix as. 



" 1 


2 ■ 




" 1 


1 ■ 




' 


1 " 












+ 






CO 


_ 




1 





2 






(30) 



The Figure [T] present the numerical errors between the exact and the numerical 
solution. 

The Figure [2] present the CPU time of the standard and the iterative splitting 
schemes. 
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Figure 1: Numerical errors of the standard Splitting scheme and the iterative 
schemes with 1, . . . , 6 iterative steps. 

Remark 3. We see that the errors decrease significantly with increasing order of 
approximation for the initialization process with the Zassenhaus splitting. Here 
we have the benefit in the application of the Zassenhaus product schemes to the 
standard Lie- Trotter or Strang-Marchuk splitting. Similar results are given with the 
iterative steps i = 3 . . . , 6, as expected. 

3.2 One phase example 

The next example is a simplified real-life problem for a multiphase transport-reaction 
equation. We deal with mobile and immobile pores in the porous media, such 
simulations are given for waste scenarios. 

We concentrate on the computational benefits of a fast computation of the mixed 
iterative scheme with the Zassenhaus formula. 
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AB, Strang, one-side 





Figure 2: CPU time of the standard Splitting scheme and the iterative schemes 
with 1, . . . , 6 iterative steps. 

The one phase equation is given as: 

dtci + V • Fci = -Aici, in n x [0,t], (31) 

dtC2 + V • Fc2 = AiCi - A2C2, in X [0, t], (32) 

F = V - DV, (33) 

ci(x, t) = ci,o(x), C2(x, t) = C2,o(x), on Q, (34) 

ci(x, t) = ci,i(x, t), C2(x, t) = C2,i(x, t), ondftx[0,t]. (35) 

where we have the parameters: v = 0.1, D = 0.01, Ai = A2 = 0.1. 

In the following we deal with the finite difference schemes for the convection and 
diffusion operators and semidiscretize the equation, which is given as: 

dtc = {Ai + A2)c, (36) 

We obtain the two matrices and consider to decouple the diffusion and convection 
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part: 



^1 



Adiff 
Adiff 



1)27x27 



Aconv 
Aconv 



+ 



-Ai 
Ai -A2 



1,27x27 



For the operator Ai and A2 we apply the sphtting method, given in Section 
The submatrices are given in the following: 



A = 



Adiff 
Adiff 

( -2 1 



^ccmv 

A 



D 





conv 
\ 



V 

Ax 



\ 

( 1 
-1 1 



V 



1 -2 1 
1 -2/ 

\ 



-1 1 

-1 1/ 



1,7x7 



where / is the number of spatial points and Ax is the spatial step size. 



/ Ai 
Ai 



Ai = 





Ai 
Ai / 



1,7x7 



Ao = 



/ A2 
A2 



V 



A2 
A2 y 



We have the following results: 

We have the spatial step size Ax = 0.1. 



d7x7 
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The Figure [3] present the numerical errors between the exact and the numerical 
solution. 

The Figure |4] present the CPU time of the standard and the iterative splitting 
schemes. 

Remark 4. For the iterative schemes with embedded Zassenhaus products, we can 
reach faster and more improved results. With 4—5 iterative steps we obtain more ac- 
curate results as we did for the expensive standard schemes. With one-side iterative 
schemes we reach the best convergence results. 

4 Conclusions and Discussions 

In this work, we have presented a novel splitting scheme combing the ideas of iter- 
ative and sequential schemes. Here the idea to decouple the expensive computation 
of only matrix exponential based schemes to simpler embedded Zassenhaus schemes, 
which have their benefits of less computational time, while the commutator can be 
computed very cheap. On the other hand simple linear iterative steps can be done 
very cheap and accelerated the solvers. The error analysis presented stable methods 
for the higher order schemes. In the applications, we could show the speedup with 
the Zassenhaus enhanced methods. In future we concentrate on linear and nonlinear 
matrix dependent scheme, that switches between no- iterative and iterative schemes 
based on Zassenhaus products. 
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